#include <assert.h>
#include <algorithm> // Required for sort()
#include "remez_exp.h"
#include "remez_log.h"
#include "remez.h"
#include "linearEquationSolver.hpp"

using namespace mpfr;

#define DEFAULT_PRECISION 224

// 若　f(x)= c0 + c1*x + c2 * x^2 + c3 * x^3 + c4 * x^4 + c5 * x^5　＋　... 是一个degree次多项式
// 多项式的系数用coeffs来表示，coeffs=[c0,c1,c2,c3,c4,c5,...], 共degree+1项
// 计算多项式在x处的值
mpreal evaluatePolynomial( const std::vector<mpreal> &coeffs, mpreal x, int degree) // degree为多项式的阶数
{
    assert(degree+1<=coeffs.size());
    mpreal r= mpreal("0");
    for (int i = 0; i <= degree; i++)
    {
        r +=  coeffs[i] * pow(x, i);
    }
    return r;
}
// 计算多项式在x处的导数值
mpreal polynomial_derivative( const std::vector<mpreal> &coeffs, 
    mpreal x, int degree)
{
    assert(degree+1<=coeffs.size());
    mpreal r= mpreal("0");
    for (int i = 1; i <= degree; i++)
    {
        r = r + mpreal(i)* coeffs[i] * pow(x, i-1);
    }
    return r;   
}

// 计算切比雪夫节点
void chebyshev_nodes(std::vector<mpreal> &nodes, mpreal a, mpreal b, int n, const char *MPREAL_FORMAT )
{
    nodes.resize(n);

    mpreal pi  = mpfr::const_pi(DEFAULT_PRECISION);
    mpreal half= mpreal("0.5");

    {
        std::string s= a.toString(MPREAL_FORMAT);
        std::cout << "a = " << s << std::endl;
    }
    {
        std::string s= b.toString(MPREAL_FORMAT);
        std::cout << "b = " << s << std::endl;
    }

    for (int i = 0; i < n; i++)
    {
        mpreal theta = mpreal(2*i+1)/mpreal(2*n) * pi;    // theta = (2.0 * i + 1) / (2.0 * n) * M_PI;
        nodes[i]= half*(a+b) + half*(b-a) * cos(theta);   // nodes[i] = 0.5 * (a + b) + 0.5 * (b - a) * cos(theta);
        {
            std::string s= nodes[i].toString(MPREAL_FORMAT);
            std::cout << "nodes[" << i << "]=" << s << std::endl;
        }

    }
}

void init_nodes_v2(std::vector<mpreal> &nodes, mpreal a, mpreal b, int n )
{
    nodes.resize(n);

    mpreal step= (b-a)/mpreal(n);
    for (int i = 0; i < n; i++)
    {
        nodes[i]=(a + step * ( mpreal(i) + mpreal(0.5) ));
    }
}

// 在区间[begin,end]之间查找极值点，使用二分查找算法,degree为多项式的阶数
mpreal search_extre_point(
    const std::vector<mpreal> &coeffs,
    Error_FP fp,
    int degree,  mpreal begin, mpreal end)
{
    mpreal low=begin;
    mpreal high=end;
    mpreal mid;

    mpreal drv_low  = fp(coeffs,low, degree); //误差函数在low处的导数值
    mpreal drv_high = fp(coeffs,high,degree); //误差函数在high处的导数值

    int count=0;
    while (low<high && count<20) // 最大迭代20次
    {
        mid = (low+high)*mpreal("0.5");
        mpreal drv_mid = fp( coeffs,mid, degree); // 误差函数在mid处的导数值

        if ( drv_mid== mpreal("0"))     //导数为0,找到极值点mid
            return mid;

        if ( drv_mid * drv_low > mpreal(0)) // mid点与low点的导数同号
        {
            low=mid;
            drv_low = fp(coeffs,low,degree);
        }
        else // mid点与high点的导数同号
        {
            high=mid;
            drv_high = fp(coeffs,high,degree);
        }
        count++;
    }
    return mid;
}

// 若fun是exp函数或者log函数,f(x)是fun(x)的近似多项式，求my_fun(x)-fun(x)在区间[a,b]之间的极值点，
// f(x)= c0 + c1*x + c2*x^2 + c3*x^3 + c4*x^4 + c5*x^5,
// degree是多项式的阶数，coeffs为多项式的系数=[c0,c1,c2,c3,c4,c5], 共degree+1个系数
// 返回极值点的, result 是极值点，i=[0..degree]，至少NPOINTS-2个极值点，至多返回至少NPOINTS个极值点
std::vector<mpreal> calc_error_extremum( \
    const std::vector<mpreal> &coeffs,\
    Error_FP fp,
    mpreal a, mpreal b, int degree,
    const char *MPREAL_FORMAT )
{
    int sample_cnt=1024;

    std::vector<mpreal> x_nodes;
    std::vector<mpreal> d_nodes;

    x_nodes.resize(sample_cnt);
    d_nodes.resize(sample_cnt);

    std::vector<mpreal> result;
    for (int i = 0; i < sample_cnt; i++)
    {
        mpreal x = (mpreal(i)/mpreal(sample_cnt))*(b-a)+a;
        x_nodes[i]= x;
        d_nodes[i] = fp( coeffs, x,degree); // 误差函数在x处的导数值,共取Ｎ个点
        // {
        //     std::string s= d_nodes[i].toString(MPREAL_FORMAT);
        //     std::cout << "d_nodes[" << i << "] = " << s << std::endl;
        // }
    }

    result.clear();
    for (int i = 0; i < sample_cnt-1; i++)
    {
        if (d_nodes[i]*d_nodes[i+1]<0) // 两个相邻点的导数异号，说明在区间[a,b]之间有极值点
        {
            mpreal x = search_extre_point(coeffs, fp,degree, x_nodes[i], x_nodes[i+1]);
            {
                //std::string s= x.toString(MPREAL_FORMAT);
                //std::cout << "x="  << s << std::endl;
            }
            result.push_back(x);
            if (result.size()== coeffs.size())
                break;
        }
    }
    return result;
}

void print_coeffs( const std::vector<mpreal> &coeffs,const char *MPREAL_FORMAT )
{
    for (size_t i = 0; i < coeffs.size(); i++)
    {
        std::string s= coeffs[i].toString(MPREAL_FORMAT);
        std::cout << "C" << i << " = " << s << std::endl;
    }
}

void print_nodes( const std::vector<mpreal> &coeffs,const char *MPREAL_FORMAT )
{
    for (size_t i = 0; i < coeffs.size(); i++)
    {
        std::string s= coeffs[i].toString(MPREAL_FORMAT);
        std::cout << "nodes[" << i << "] = " << s << std::endl;
    }
}

// 从多个节点中挑选 n 个节点，使得多项式在节点上的误差最小
 std::vector<mpreal> select_nodes(
    const std::vector<mpreal> &coeffs,
    Error_FP fp1,
    Error_FP fp2,
    mpreal a, mpreal b, int degree, int n,
    const char *MPREAL_FORMAT)
{
    std::vector<mpreal>  tmp_nodes;
    std::vector<bool>  nodes_is_valid;

    std::vector<mpreal> nodes = calc_error_extremum( coeffs,fp1, a,  b, degree,MPREAL_FORMAT);

    if ( nodes.size() +2  < n) // point_cnt个节点，在加上额外的两个节点a,b, 总数不能小于NPOINTS
    {
        std::cerr << "Error, cannot get " <<  n-2 << " nodes" << std::endl;
        nodes.clear();
        return nodes;
    }

    int tmp_nodes_cnt = nodes.size()+2;
    tmp_nodes.resize(tmp_nodes_cnt);
    nodes_is_valid.resize(tmp_nodes_cnt);
    for (int i=0;i<tmp_nodes_cnt;i++)
    {
        if (i==0)
            tmp_nodes[i]=a;
        else if (i == tmp_nodes_cnt-1)
            tmp_nodes[i]=b;
        else
            tmp_nodes[i]= nodes[i-1];
        nodes_is_valid[i]=true; // mark all node are valid
    }

    std::sort(tmp_nodes.begin(), tmp_nodes.end());
    int   valid_nodes_cnt= tmp_nodes_cnt;
    while ( valid_nodes_cnt > n )
    {
        mpreal min = const_infinity(1); // positive infinity
        int idx_on_min=0;
        for (int i=0;i<tmp_nodes_cnt;i++)
        {
            mpreal abs_error = abs( fp2( coeffs, tmp_nodes[i],degree));
            
            if ( abs_error<min &&  nodes_is_valid[idx_on_min])
            {
                min=abs_error;
                idx_on_min=i;
            }
        }
        nodes_is_valid[idx_on_min]= false; // mark node as removed
        valid_nodes_cnt--;
    }

    nodes.clear();
    for (int i=0;i<tmp_nodes_cnt;i++)
    {
        if (nodes_is_valid[i])
        {
            nodes.push_back(tmp_nodes[i]);
        }
    }
    return nodes;
}

std::string get_ordinal_number( int i)
{
    switch (i)
    {
        case 1: return std::string("first"); break;
        case 2: return std::string("second"); break;
        case 3: return std::string("third"); break;
        default: return std::to_string(i) +"th";
    }
}

std::vector<double> get_glibc_coeffs()
{
    std::vector<double> coeffs;
    coeffs.resize(6);

    coeffs[0]=1.0;
    coeffs[1]=1.0;
    coeffs[2]=0x1.ffffffffffdbdp-2;
    coeffs[3]=0x1.555555555543cp-3;
    coeffs[4]=0x1.55555cf172b91p-5;
    coeffs[5]=0x1.1111167a4d017p-7;

    for (int i=1;i<=5;i++)
        printf("glibc C%d=%.18lf\n",i,coeffs[i]);

    return coeffs;
}




